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Abstract 

We apply the wavelet transform modulus maxima (WTMM) method § to the anal- 
ysis of simulated MBE - grown surfaces. In contrast to the structure function approach 
commonly used in the literature, this new method permits an investigation of the com- 
plete singularity spectrum. We focus on a kinetic Monte-Carlo model with Arrhenius 
dynamics, which in particular takes into consideration the process of thermally activated 
desorption of particles. We find a wide spectrum of Holder exponents, which reflects the 
multiaffine surface morphology. Although our choice of parameters yields small desorp- 
tion rates (< 3%), we observe a dramatic change in the singularity spectrum, which is 
shifted towards smaller Holder exponents. Our results offer a mathematical foundation of 
anomalous scaling: We identify the global exponent a g with the Holder exponent which 
maximizes the singularity spectrum. 



1 Introduction 

Inspired by the great technological importance of epitaxial crystal growth, the past decade has 
raised much theoretical research in the subject of kinetic roughening of surfaces during growth. 
The investigation of this effect, which is undesirable in practical applications, promises deep 
insight into statistical physics far from thermal equilibrium, see e.g. 0| for an overview. 
We focus on a full-diffusion Monte-Carlo model of homoepitaxial growth of a hypothetical 
material with simple cubic lattice structure under solid on solid conditions, i.e. the effects of 
overhangs and displacements are being neglected. Then, the crystal can be described by a 
two-dimensional array of integers which denote the height f(x) of the surface. On each site, 
new particles are deposited with a rate r a . Particles on the surface are hopping to nearest 
neighbour sites with Arrhenius rates ^o ex P( — (Eb + nE n ) / \k\,T)) , where E b and E n are the 
binding energies of a particle to the substrate and to its n nearest neighbours, uq is the attempt 
frequency, and k b T has its usual meaning. In contrast to earlier investigations of similar 
models ||, we permit the desorption of particles from the surface with rates vq exp(— (Ed + 
nE n )/(k b T)), where E d > E b . 

The aim of this publication is twofold: We will first discuss the advantages of the wavelet 
analysis compared to the structure function (SF) approach, which has to date solely been used 
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in the investigation of multiaffine surfaces. Then, we will apply this formalism to investigate 
the influence of desorption on kinetic roughening. We conclude with some remarks on the 
relevance of universality classes for our results. 



2 Scaling concepts 

The standard approach of dynamic scaling |Q| assumes that the statistical properties of a 
growing surface before saturation remain invariant under a simultanous transformation of 
spatial extension x, heigth f(x) and time t, 

£^bx = x; f^b a f = f; t^b z t = t'; (1) 

where b is an arbitrary positive constant. This implies, that a part of the surface smaller 
than the correlation length ~ t l l z can be regarded as self-affine with Hurst exponent a. 
A popular method of measuring a uses heigth-heigth correlation functions of (theoretically) 
arbitrary order q: 

G(q, I, t) := (\f(x, t) - f(x + I, t)\ Q ) £ ~ l qa g(l/m), (2) 

where g(x) — > const, for x — > and g{x) — > const. • x~ qa for x — > oo. In practice, q = 2 is the 
most common choice. 

In principle, there are two different ways to measure a: The local approach determines a 
from the initial slope of ln(G(q, I, t)) versus ln(Z) for small I. The global approach analyzes 



the dependence of the surface width w = y((f(x,t) — {f(x,t))) 2 )g in the saturation regime 
on the system size N: w sa t(N) ~ N ag . Before saturation, the surface width increases like 
w ~ t@ \ where (3 = a/z. An alternative which avoids the simulation of different system sizes 
uses the complete functional dependence of equation g a g and z are chosen such that the 
curves of G(2,l,t)/l 9 versus ljt x l z collapse on a unique function g within a large range of t 
and /. 

However, a careful analysis of simulation data [Q, ||, ||, |3| has shown, that several models 
of epitaxial growth show significant deviations from this simple picture. First, one obtains 
different values of a from the local than from the global approach, a phenomenon which 
is called anomalous scaling. Second, one often finds multiscaling: height-height correlation 
functions of different order yield a hirarchy of (/-dependent exponents a(q), when determined 
from the initial power-law behaviour of G(q, I, t). 

These observations can be interpreted within the mathematical framework of multifrac- 
tality: The Holder exponent |l|, [2L [l^, |l(J h(x~o) of a function / at xq is defined as the largest 
exponent such that there exists a polynomial of order n < h(xo) and a constant C which 
yield \ f(x) — P n (x — xq)\ < C\x — xo\ h ( x °^ in the neighbourhood of xq. The Holder exponent 
is a local counterpart of the Hurst exponent: a self-affine function with Hurst exponent a 
has h(x) = a everywhere. However, in the case of a multiaffine function different points x 
might be characterized by different Holder exponents. This general case is characterized by 
the singularity spectrum D(h), which denotes the Hausdorff dimension of the set of points, 
where h is the Holder exponent of /. 
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3 The wavelet approach to multifractality 



There is a deep analogy between multifractality and thermodynamics [||, |?], ||, where the 
scaling exponents play the role of energy, the singularity spectrum corresponds to entropy, 
and q plays the role of inverse temperature. So, theoretically D{h) might be calculated via 
a Legendre transform of a(q): D(h) = mm q (qh - qa(q) + 2) Q, fE^, 0, a method which 
has been called structure function (SF) approach. However, its practical application raises 
fundamental difficulties: First, to obtain the complete singularity spectrum, one needs ct(q) 
for positive and negative q. But as \f(x,t) — f(x + l,t)\ might become zero, G(q,x,t) is in 
principle undefined for q < 0. Therefore, only the left, ascending part of D(h) is accessible to 
this method. Additionally, the results of the SF method can easily be corrupted by polynomial 
trends in f(x) |I2[ . It might be due to these difficulties, that - to our knowledge - no attempt 
to determine the singularity spectrum of growing surfaces from a(q) has ever been made. 
Although it has been argued that the a(q) collapse onto a single a in the limit t — > oo, which 
characterizes the asymptotic universality class of the model |3], |6|, |ll]] , we are convinced that 
deeper insight into fractal growth on experimentally relevant finite timescales can be gained 
from a detailed knowledge of the D(h) spectrum. 

To this end, we follow the strategy suggested by Arneodo et. al |Q, ||, 12], which circum- 



vents the problems of the SF approach and permits a reliable measurement of the complete 
D(h). Mathematically, the wavelet transform of a function f{x) of two variables is defined as 
its convolution with the complex conjugate of the wavelet ip, which is dilated with the scale 
a and rotated by an angle 9 ||: 

T^[f](b,6,a) = C- 1/2 a- 2 J d 2 x ^(a^R-^x - b))f(x). (3) 

Here R# is the usual 2-dimensional rotation matrix, and CU = (2ir) 2 J d 2 k\k\~ 2 \ip(k)\ 2 is a 
normalization constant, whose existence requires square integrability of the wavelet ip{x) in 
fourier space. Apart from this constraint, the wavelet can (in principle) be an arbitrary 
complex- valued function. Introducing the wavelet ips(x) = 5(x) — 5(x + n), where n is an 
arbitrary unit vector, one obtains easily 

T^ s [f](jb,e,a) =C~ l J 2 [f(b)-f(b + aRgn)\ => J d 2 b \t^ 5 [f}(b, 6, a)\ 9 oc G(q, align). (4) 

Consequently, a calculation of the moments of the wavelet transform of the surface yields 
the SF approach as a special case. To avoid its weaknesses, two major improvements are 
necessary: 

First, we use a class of wavelets with a greater number of vanishing moments than 
i[)${x). This increases the range of accessible Holder exponents and improves the insensitivity 
to polynomial trends in f(x). We introduce a two-component version of the wavelet transform 

where the analyzing wavelets ^ are defined as partial derivatives of a radially symmetrical 
convolution function <I>(x): ^i(x) = d^/dx, ^(x) = d<&/dy. Then T^[f](b, a) can be 
written as the gradient of f(x), smoothed with a filter <3? w.r.t. b. This definition becomes 
a special case of equation [|, when multiplied with fig = (cos(#), sin(#)) T , yet allows for 
an easier numerical computationF]. For example, $ can be a gaussian, where = 1, or 
For simplicity, the irrelevant constant has been omitted 
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3>i(x) = (2 — x 2 ) exp(— x 2 /2), which has two vanishing moments. 

Second, the integration over b in equation || is undefined for q < 0, since the wavelet 
coefficients might become zero. The basic idea is to replace it with a discrete summation 
over an appropriate partition of the wavelet transform which obtains nonzero values only, 
but preserves the relevant information on the Holder regularity of f{x). In the following, we 
will give a brief outline of the rather involved algorithm and refer the reader to [jl], |9], [l2j for 
more details and a mathematical proof. The wavelet transform modulus maxima (WTMM) 
are defined as local maxima of the modulus M$[f](b, a) := [/](&, a)\ in the direction 
of T^[f](b, a) for fixed a. These WTMM lie on connected curves, which trace structures 
of size ~ a on the surface. The strength of each is characterized by the maximal value 
of M^[f](b, a) along the line, the so-called wavelet transform modulus maxima maximum 
(WTMMM) [|]]. While proceeding from large to small a, successively smaller structures are 
resolved. Connecting the WTMMM at different scales yields the set C of maxima lines /, 
which lead to the locations of the singularities of f(x) in the limit a — > 0. The partition 
functions 

Z(q,a)= J2 ( SU P M f[/](M0] ~« r(9) fora^0 (6) 

are defined on the subset C(a) of lines which cross the scale a. Prom the analogy between 
the multifractal formalism and thermodynamics, D{h) is calculated via a legendre transform 
of the exponents r(q), which characterize the scaling behaviour of Z(q,a) on small scales a: 
D(h) = mm q (qh — r(q)). Additionally, r{q) itself has a physical meaning for some q: — r(0) 
is the fractal dimension of the set of points where h{x) < oo, while the fractal dimension of 
the surface f(x) itself equals max(2, 1 — r(l)). 



4 Results 

In our simulations, we choose the parameters uq = 10 12 /s, E\, = 0.9 eV and E n = 0.25 eV, and 
a temperature T = 450K. To study the influence of desorption, we consider three models with 
different activation energies Ed- in model A desorption is forbidden, i.e. Ed = oo. Models B 
and C have = 1.1 eV and Ed = 1.0 eV. We simulate the deposition of 2 • 10 4 monolayers 
at a growth rate of one monolayer per second on a lattice of N x N unit cells using periodic 
boundary conditions, our standard value being iV = 512. To check for finite size effects, we 
have also simulated iV = 256. In all presented results averages over 6 independent simulation 
runs have been performed. Although we have used an optimized algorithm, these simulations 
consumed several weeks of CPU time on our workstation cluster. 

First, we have checked our results for artifacts resulting from properties of the analyzing 
wavelet rather than from the analyzed surface by using different convolution functions 
$o is the gaussian function, $ n ,n > 1 are products of gaussians and polynomials, which have 
been chosen in a way that the first n moments vanish. Then, the analyzing wavelets have 
n $ =n + l vanishing moments. We find (figure [l] a), that the r(g)-curve obtained with 

deviates significantly from those obtained with $i, $2 and $3. The latter agree apart 
from small differences which are mainly due to the discrete sampling of the wavelet in the 
numerical implementation of the algorithm. This is explained by the theoretical result [|Toj| 
that dr(q)/dq = for q < q cr u. < if the number of vanishing moments of the analyzing 
wavelet is too small. Consequently, the agreement of the other curves proves their physical 
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Figure 1: Panel a: r(g) in the = oo - model as obtained from investigations with different 
wavelets. Note the deviations of the data obtained with the gaussian function <I>o- Panel 
b: r{q) for models with different activation energies of desorption. All data have been 
obtained from surfaces after 2 • 10 s of simulated time. Sizes of errorbars are on the order of 
symbol sizes. 

relevance. 

Figure |] b shows averages of r(q) curves obtained with the convolution functions $i, 
q?2, $3 from surfaces after 2 • 10 4 s of growth on an initially flat substrate. For all our 
models, their nonlinear behaviour reflects the multiafhne surface morphology. From the fact 
that these curves are reproduced within statistical errors in simulations with N = 256, we 
conclude that finite size effects can be neglected. Clearly, desorption reduces the slope of r(q), 
although only a small fraction of the incoming particles is desorbed: 0.18% in model B and 
2.57% in model C with slightly higher values at earlier times. The corresponding singularity 
spectra are shown in figure ^ a. They have a typical shape whose descending part seems to 
be symmetrical to the ascending part and which changes at most slightly, while the whole 
spectra are shifted towards smaller Holder exponents as desorption becomes more important. 
We emphasize that we find no evidence for a time dependence of the singularity spectra within 
the ranee 9700s < t < 2 • 10 4 s, so that our results do not support the idea of an asymptotic 
regime characterized by a single exponent a. However, the accessible time range of computer 
simulations is limited, so that we cannot finally disprove the existence of such a regime. 

The multifractal formalism has replaced the unique scaling exponent a of spatial exten- 
sion in the simple picture of dynamic scaling (equation |l|) with a wide spectrum of Holder 
exponents. By analogy, one might find it necessary to replace the scaling exponent (3 with 
a distribution of temporal counterparts of h. To answer this question, we investigate the 
propability distribution function (PDF) P(f — (/) , t) of surface heights. Dynamical scale 
invariance with a single (3 demands that 

P(f-(f},t) = P( L ^ 11 )^ (7) 

i.e. the rescaled PDFs Pt@ should collapse onto a single function P when plotted as a function 
of (/ — (f))/t@ within a large time range. 




Figure 2: Panel a: singularity spectra obtained from a Legendre transform of the data in fig. 
H b. Panel b: Data collaps of rescaled PDFs of surface heights for times between 150s and 
2 • 10 4 s (model A) respectively 150s and 7500s (model C). We have used /3 = 0.188 for model 
A with Ed = oo and (3 = 0.109 for model C (Ed = l.OeV). Time is measured in seconds. 

We measure (3 from the increase of the surface width with time, which follows a power 
law for t > 150s in models A and B ((3 A = 0.19 ± 0.01, (3 B = 0.17 ± 0.01) respectively 
150s < t < 7500s in model C ((3c = 0.11 ± 0.01), which then starts to approach the final 
saturation regime. The high quality of the data collapse of the PDFs shown in figure ^ b 
proves that the scaling form [?] holds, showing that a single exponent describes the scaling 
behaviour of P(f — (/) , t). This parallels the finding of Krug in || for the one-dimensional 
Das Sarma-Tamborenea model. 

Finally, the WTMM method, which is a precise tool to investigate local scaling properties 
of surfaces might help to get some insight into the phenomenon of anomalous scaling. The 
conventional picture []l3| , O, 15 1 notes the difference between the global a g and a "local a" 



which is determined from the power-law behaviour of G(2,l,t) for small and, within the 
multifractal formalism, simply corresponds to a Holder exponent on the ascending part of 
the singularity spectrum. We have determined the global scaling exponents a g and z from 
the data collapse of the scaled height-height correlation function G(2, 1, t) and find agreement 
within statistical errors between a g and that value of the Holder exponent h m which maxmizes 
D(h) (table |]). This empirical result can be explained with a saddle-point argument: We 
calculate the surface width 

w2 = w*j d2x{f{s) " (/>)2 = w^ dl ! d2x 6{ ~ h " - » 2 • ( 8 ) 



1(h) 



Since 1(h) grows like N D ^ with the system size, in large systems the integral over h will 
be dominated by I(h m ). That means, that w and therefore the global scaling properties of 
the surface are governed by the subset of points, which has the greatest fractal dimension. 
Consequently, the surface will behave like a self-affine surface with Hurst exponent h m on 
lenghtscales comparable to the system size. 
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Model 


Pd 







a g 


Z 9 


Df 


A 





0.19 ± 0.01 


0.54 ± 0.01 


0.55 


2.9 


2.32 ± 0.01 


B 


0.18 % 


0.17 ± 0.01 


0.52 ± 0.01 


0.51 


3.3 


2.35 ± 0.01 


C 


2.57 % 


0.11 ± 0.01 


0.38 ± 0.01 


0.39 


3.5 


2.45 ± 0.02 



Table 1: Simulation results: pd is the fraction of particles which desorbes, f3 is the scaling 
exponent of the PDF of surface heights, h m the Holder exponent which maximizes D(h), a g 
and z g are the global scaling exponents of (7(2, l,t), Df is the fractal dimension of the surface. 

5 Conclusions 

Table [l] summarizes our results. Model A without desorption reviews the results in [||, which 
have been obtained with slightly different activation energies on smaller systems and shorter 
timescales. Models B and C show, that desorption is an important process, which, although 
it affects only a small fraction of the adsorbed particles, must not be neglected, since it alters 
the scaling properties of the surfaces by reducing j3 and by shifting the singularity spectrum 
towards smaller Holder exponents. Since the scaling behaviour depends strongly on the height 
of the energy barrier of desorption, and the singularity spectra have no measurable tendency 
to narrow with time, our results can not be used to make any decision on the aymptotic 
universality class of the investigated model. However, they show, that the paradigm of a 
few universality classes characterized by a small number of exponents, which are indepen- 
dent on details of the model, is not adequate to catch the features of kinetic roughening on 
experimentally relevant timescales of a few hours of growth. 

We are convinced, that the application of new mathematical tools like the wavelet analysis 
will help to find a better description of fractal growth phenomena in the future. 

We thank A. Arneodo and J. M. Lopez for providing us recent preprints before publication 
and A. Freking for a critical reading of the manuscript. 
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